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Abstract. Intense-field ionization of the hydrogen molecular ion by linearly-polarized light 
is modelled by direct solution of the fixed-nuclei time-dependent Schrodinger equation and 
compared with recent experiments. Parallel transitions are calculated using algorithms which 
exploit massively parallel computers. We identify and calculate dynamic tunnelhng ionization 
resonances that depend on laser wavelength and intensity, and molecular bond length. Results 
for A ~ 1064 nm are consistent with static tunnelling ionization. At shorter wavelengths 
A ~ 790 nm large dynamic corrections are observed. The results agree very well with recent 
experimental measurements of the ion spectra. Our results reproduce the single peak resonance 
and provide accurate ionization rate estimates at high intensities. At lower intensities our 
results confirm a double peak in the ionization rate as the bond length varies. 



The mechanism of high-intensity ionization by infrared and optical wavelength light is often 
considered a static tunnelling process. The simplicity of this model is hugely appealing 
because of the ease of calculation. The ionization rates are effectively independent of 
wavelength, and to some extent the internal structure of the molecule can be ignored 1 1 1. A 
rough criterion for validity of this model is given by the Keldysh parameter, 7^ = yJ\Ei\/2Up, 
where the internal binding energy is {\Ei\) and the external laser-driven kinetic energy is 
{Up). When the conditions are such that 7^ <IC 1, the ionization process is dominated 
by static tunnelling in which the shape of the potential strongly (exponentially) affects the 
ionization rate. At certain critical distances between the nuclei, discovered by Codling and 
co-workers |2J, the ionization rate can rise sharply producing a sequence of fast fragments 
ions at sharply-defined energies. Predictions for ion yields and energies based on classical 
arguments f\\ agree very well with experiments even for large diatomic molecules such as I2. 
The presence of critical distances would be evident in polyatomic molecules and is also seen 
in small rare-gas clusters |3|. The tunnelling process is generally relatively fast compared 
to the vibrational motion of the molecule, so the fixed-nuclei approximation is reasonable. 
However the tunnelling time may be longer than the optical period of the laser Under these 
conditions the process is more accurately termed a dynamic tunnelling process. In this paper 
we provide evidence of just such effects for Ti:Sapphire light A ^ 790 nm at intensities 
I ~ 10^^ W cm^^. Our theoretical results in this wavelength region do not agree with cycle- 
averaged static field models, however the results do agree well with the features observed in 
experimental studies. 

In a molecule with few electrons the ionization process can be studied quantally with few 
approximations. For one-electron models, static-field ionization resonances in the potential 



Letter to the Editor 



2 



wells |4l |5J |6l occur at distances far from the equilibrium internuclear separation and 
tend to produce low-energy ions. Experiments have confirmed the existence of enhanced 
multiphoton ionization in the hydrogen molecular ion at infrared wavelengths |9, 10 1 but 
at intensities such that dynamic effects of the field cannot be neglected | IJJ . The well- 
established Fourier-Floquet analysis (6| 1121 is not particularly suitable for the study of 
long-wavelength excitations as the number of frequency components required is very large. 
Moreover, this approach supposes continuous wave conditions such that the state decays 
exponentially from an isolated resonance state with a lifetime longer than the optical cycle 
or natural orbital period. Conversely, long-wavelength pulses can be described by quasistatic 
fields under the conditions such that 7^ ^ 1. However, simple tunnelling formulae assume 
exponential decay from a single isolated resonance connected adiabatically to the field-free 
state. This neglects nonadiabatic transitions within the well |7, 13 1 and rescattering of the 
continuum electron I14II15I . Given these difficulties, the direct solution of the Schrodinger 
equation has distinct advantages. It is suited for all intensities and electronic states and all 
wavelengths and pulse shapes. In particular it is capable of describing pulse-shape effects and 
nonadiabatic transitions. Thus it is highly appropriate for realistic modelling of experiments 
at infrared wavelengths such that 7^ ~ 1. 

Tunable Ti:Sapphire light A ^ 780 — 800 nm interacting with atomic hydrogen for 
example, achieves the pure tunnelling regime, 7^ « 0.1, only for intensities / > 1 x IG^^W 
cm-2, while for A = 800 nm with / 3 x lO^-^W cm"^ |9|^ ^ g.T. Under these latter 
conditions the ionization rate is well-defined, however a static tunnelling model is unlikely 
to give a correct estimate of the resonance positions and rates; we show that this is indeed 
the case. In fact, the dynamic effects displace the critical distances, change the ionization 
rates, and create electron excitation resonances. In the present work we solve the electronic 
dynamics exactly by a direct solution of the time-dependent Schrodinger equation (TDSE). 
This does not include broadening due to the finite focal volume and the corrections due to 
nuclear motion. Using atomic units, the ground state of the molecular ion is characterized 
by a bond length = 2.0 and rotational constant Be = 1.36 x 10"'*. If the laser pulse 
duration is relatively short 20 fs compared with rotational timescale for this molecule 
{Trot ~ l/^o ~ 200 fs), then the laser-molecule interaction can be regarded as sudden in 
comparison to the rotation of the system. Neglect of rotation effects is therefore reasonable. In 
spite of these simplifications the results are very promising and in remarkably good agreement 
with experiment for the ion energy spectrum. The indications are that appropriate refinements 
of the model would improve agreement, but this remains a goal for future work. 

Making these approximations the TDSE, in atomic units, reads 

iJe*.(r-,t) = i^^',(r,t), (1) 
ot 

where is the electronic Hamiltonian and (r,<) the electronic wavefunction. 
Monochromatic light with linear polarization parallel to the internuclear axis implies a 
cylindrical symmetry about the internuclear axis with the associated good quantum number 
A. Thus the electron position can be completely described by the radial, p, and axial, z, 
coordinates with respect to an origin taken at the midpoint between the nuclei. The TDSE 
reduces to a 2h-1 dimensional partial differential equation where the electronic wavefunction 
is written as '^e{p, z, t) and the electronic Hamiltonian takes the form 

H,{R;p,z;t)^~- + _ + -_ j + _ + p, z) + t), (2) 

where R is the distance between the two nuclei which have charges Zi and Z2, Vc{R, p, z) is 
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Internuclear separation, /?(a.u.) 



Figure 1. Ionization rate, F, for dynamic and cycled-averaged-static fields. Tlie figure shows 
the wavelength and bond length dependence of the ionization rate. The maximum field strength 
Eq = 0.05338 a.u. corresponds to / = 1x10^"' W cm . Time-independent cycle- 
averaged static field rates |6|, time-dependent rates for A = 1064 nm (this work), O ; 
time-dependent calculations A=790 nm (this work). A; time-dependent rates for A = 1064 
nm m, • . 



the electronic potential given by 

Vc[R,p,z) = - 



(3) 



and Vn\-i{z, t) the molecule-laser interaction in the length gauge is given by 

Vn,^i{z, t) = zEQf{t) COS wi, (4) 

where Eq is the peak electric field, lo the angular frequency and f{t) is the pulse envelope 
given by 



1 



fit) = 



(S) 



(5) 



< i < Ti 

Tl <t <,Tl+T2 

t < 0, t > r2 + 2ti 

where the pulse ramp time is ti and the pulse duration T2, with associated bandwidth 
Aw — 1/t2. In the calculations presented in this paper ri = 5 cycles and T2 = 10 cycles. It is 
convenient to change the dependent variable to remove the first-derivative in p as follows 1161 

<i>(p,z,t)^{2npY'H,{p,z,t), (6) 

so that for E-symmetry (A = 0) the time-dependent equation is 

[T, +Tp + Kn(p, z, R) + Kn_i(z, t)] 4>{p, z, t) = i^0(p, z, t) , (7) 



where 



1 



dz"^ 



(8) 
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This 2+1 dimensional TDSE can be discretized on an Np x Nz x Nt space-time grid. We 
label the Np radial grid points by, {pi, p2, . . . pi, . . . pn }, while the axial grid points 
are denoted by, {zi, Z2, ■ ■ ■ Zj, . . . z^^}. The time evolution progresses through the sequence 
of times {ti, t2, ■ ■ - tk, ■ ■ ■ tjVt}- In this case the wavefunction can be written as the array 
4>{p, z, t) (t){zi, pj,th). The method of discretization of the Hamiltonian divides the axial 
and radial coordinates into subspaces. Two distinct but complementary grid methods are used 
for the subspaces 1 1 6 1 . The radial subspace is discretized on a semi-infinite range using a small 
number Np of unevenly spaced points that are the nodes of global interpolating functions; 
Lagrange meshes. This leads to a small dense matrix for the Hamiltonian in the p-subspace. 
On the other hand the axial coordinate subspace is represented by a large number of equally- 
spaced points, with spacing Az = 0.1 a.u., as lattice points of a finite-difference scheme. 
The associated subspace Hamiltonian matrix is large but sparse. Our approach is tailored 
to the requirements of accuracy and computational efficiency. This approach can easily be 
parallelized to make use of massively parallel processors 1141 1161 . At the very least, the 
dimensions of the cylindrical box, height 2z,„ax radius pmax, must be chosen to encompass 
the tightly-bound states of the system. At the same time the box should be large enough to 
allow the continuum states to evolve unfettered. As the wavefunction approaches the edge of 
the box boundaries, we capture the photoelectrons by employing a masking function to absorb 
the outgoing flux L17J . The ground state is calculated via an iterative Lanzcos calculation 
as described in 1161 . 

The quasistatic nature of long-wavelength pulses (A ~ 1064 nm) means that it is fair to 
compare the cycle-average static field ionization rate with the time-dependent ionization rate 
fTl . The dynamic-field (wavelength-dependent) effects can be judged from figureQin which 
we choose the wavelengths A — 79Gnm and A = lG64nm with the same average intensity 
I = 1 X 10^** W cm~2 (Eq = 0.05338 a.u.). Firstly, for A = 1064nm the cycle-averaged static 
field features 0121 are very similar to those found using our time-dependent method; with 
two resonance peaks near R ^ 5 a.u. and i? ~ 9 a.u. However, there are large differences 
in the shape and relative heights of the peaks. The prominent resonance near R ^ 9 a.u. is 
the charge-resonance peak EJ. The inner peak (i? ^ 5 a.u.) is a feature of the potential 
barrier. The longer wavelength A — 1064 nm does give results very similar to the static cycle- 
averaged results as expected [6 1 . For A ~ 790 nm, there is a significant reduction in the heights 
of these peaks and some indication of peak positions moving towards smaller internuclear 
separations. Calculation at 390nm demonstrate that this trend in peak position moving to 
smaller values of R is maintained with the first peak found at i? = 4 a.u. for this wavelength. 
The resonance structure depends on both bond length and wavelength. It is interesting that as 
the molecule separates into its atomic fragments, the wavelength dependence disappears and 
the static field result is valid. Figure [2 illustrates very clearly that molecular field ionization 
differs strongly from atomic field ionization for the bond lengths, wavelengths and intensities 
of interest. Indeed the molecular ionization rates only converge to within 5% of the atomic 
rates at i? = 20 a.u. For instance at a wavelength of 1064nm the molecular ionization rate 
at i? = 20 is 2.85 x 10"^ fs^^ compared with the atomic rate of 2.92 x 10"^ fs^^. These 
atomic rates are calculated using the present code which takes R=Q and Z\ = ^ 
These atomic results are in agreement with other accurate time-dependent results |8 1 to within 
0.5%. Our time-dependent results shown in figure^Jfor A = 1064 nm are consistent with 
previous static field cycle-averaged results |6 7|, although these results disagree with other 
time-dependent results L5J. There are strong similarities in the R dependence of the rates 
with those calculated previously. However our rates are up to 4 times higher than those of (Sj 
and the resonance positions are displaced to smaller R values, so that we predict faster ion 
fragments with higher yields. 
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Internuclear separation, R (a.u.) 

Figure 2. Comparison of ionization rates with experimental measurements at A=800 nm, 
cm ^. Tlie experimental results |9| (• ) have been normalized to the 
theoretical calculations (A). 

We apply our model to simulate experiments (9l on the ion energy spectra from the 
dissociative ionization of the hydrogen molecular ion using A ~ 800 nm at / = 3.2 x 10^^ W 
cm^^; the results are presented in figure|2] In the experiment an H2 gas was ionized to form 
the ion target. The molecular ions subsequently dissociate and ionize in the laser field, 
with the proton fragments extracted and energy analyzed. By relating the kinetic energies 
of the ions to the Coulomb explosion curve, ionization rates could be deduced for the range 
of molecular bond lengths. The experimental conditions were such that saturation of the 
ionization channel was avoided, permitting ion yields from larger R values to be estimated 
and eliminating ion yields from the larger focal volume. Our approximation in assuming a 
very localized high intensity region, smaller than the diffraction limit, is well justified. The 
sensitivity of the ionization rate to changes in intensity and wavelength was noted in the results 
presented in figure ^ Consider the changes in the ionization rates in going from the data 
presented in figure^for A = 790 nm to that for A = 800 nm at an intensity 3.2 times higher, 
in figure |2] The ionization rates in figure |2] are roughly 15 times higher which is consistent 
with an exponential increase giving rise to more easily measurable ion yields. However the 
double peak structure of figure[l]is now dominated by a single broad maximum near i? ^ 7. In 
comparing with experiment we have in figure|2lnormalized the laboratory data to our results. 
We see that the shape of the theory and experimental curves are in remarkable agreement, 
in spite of the assumptions made. The single broad peak is reproduced rather well, although 
some additional structure present in the simulations is not resolved by experiment. For R > 8 
the theory and experiment are in disagreement, and the theoretical estimate of ion yield from 
i? ^ 9 is much lower than the experimental results. This might be attributed in part to the 
variation in focal volume intensity. At the edges of the focal spot the intensities decrease 
but the interaction volume is larger |1|. Moreover at lower intensities the field ionization 
rates move to larger i? (6|. So one would expect that an inclusion of focal volume variation 
would broaden the peak to larger values of R and partially compensate for this shortfall. 
The second feature is that the theoretical results predict high ionization far in excess of that 
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Figure 3. Comparison of ionization rates from the present time-dependent calculations (A) at 
A=790 nm and / = 0.6 X 10^"'W cm"^ (Eq = 0.041 a.u.) with static field cycle-average 
rates \6>i (♦) at a comparable intensity of / = 0.56 X 10^"'W cm"^ Eq = 0.040 a.u. 

found experimentally at large R. The i? — > oo atomic limit of the theory results is very 
accurately known and consistent with the theoretical results in figure |2| The theoretical data 
for ionization rates can be considered as accurate. However the lower ion yield observed 
can be explained by the fact that during molecular dissociation the molecular ion can ionize 
at smaller R values. If the ionization rates are large at small R, then few if any molecules 
can survive to be ionized at large bond lengths. A rough estimate of the survival probability 
P{R) of the molecular ion can be found from the classical dynamics of the ionized molecular 
ion. The depletion rate is given by dP/dR « —{T/v)P where v{R) is the classical relative 
velocity of the protons such that uipv'^ /A + Z1Z2/R = ZiZ2/Re and nip is the proton mass. 
For the case A — 800 nm and / 3.2 x lO^"' Wcm^^ a rough estimate based on this model 
gives P « 0.2 at i? = 14. This is only an indication of the effect but it is consistent with the 
findings of Dundas 1 1 5 1 . 

Very recently data have become available from experiments in Garching on for 
A = 790nm at intensities just above the Coulomb explosion threshold, namely / = 0.6 x lO^"* 
W cm^^ 1 18|. In these results the first observations of vibrationally resolved structure has 
been obtained. From the ion momentum distribution, it was suggested that a critical distance 
around i? = 12 could explain the results. Existing static field cycle-average rates at a 
comparable intensity / = 0.56 x lO^^W cm^^, iJp = 0.04 a.u. f6l are shown in figure|3] 
These results confirm that the ionization rates are extremely small, the reduction in intensity 
by a factor of 5 leading to ion yields roughly one hundred times smaller The double peak 
structure emerges in our calculations, and since the rates are now reduced, the bulk of the 
molecules will reach the outer resonance position. In figure|3]the time-dependent calculations 
for A=790 nm and / = 0.6 x lO^'^W cm^^ are in fairly good agreement with the static field 
results. However we note (figure|3ji the inner peak i? 6 — 8 is broad and high and ought 
to produce ion yields comparable to the sharp outer peak near R — 11. The observation 
of quantal vibrational structure in the ion spectrum means that a full quantum treatment of 
nuclear dynamics is required to analyze these new experiments in full. 

We have solved the full-dimensional TDSE for the electron dynamics of HJ in linearly 
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polarized laser fields, assuming that the nuclei are fixed in space. The method employed 
is highly accurate and can be efficiently implemented on parallel processing computers. 
Ionization rates can be calculated for all nuclear separations and for wavelengths from the 
infrared to x-ray, for a range of laser pulses. Comparison of our results with other theoretical 
calculations and recent experimental measurements show very good agreement. We have 
been able to identify and calculate dynamic tunnelling resonances for A = 790 nm and 
A = 800 nm and obtain accurate estimates of the ionization rates and large dynamic tunnelling 
corrections are observed. A major simplification in the model is the fixed-nuclei assumption. 
However, our results for A = 800 nm and / = 3.2 x lO^"* W cm^^ reproduce the measured 
dependence of ionization rate on bond length. At shorter wavelength and lower intensities, 
A — 790 nm and / = 0.6 x lO'^'* W cm~^, our results indicate a double peak structure 
in the ionization rate as the bond length varies. The outer ionization resonance agrees with 
experimental measurements [T8I| . 

To model the experiments more realistically several extensions to the current approach 
are required. Firstly, the energy and angular momentum exchanges between the nuclei and 
electrons will occur during process. Dundas 1 15 1 has combined the full electronic dynamics 
with a quantal vibrational motion for intense field dissociative ionization and found that the 
dynamic tunnelling resonances dominate strongly over pure dissociation at high intensities. 
A classical model of nuclear motion will not be sufficient as the wavepacket will disperse 
during the process and indeed experiments are now able to resolve the vibrational structure in 
the ion yield 1 18 1. The quantal motion is essential to obtain an ion spectrum distribution rather 
than one-to-one mapping of ion energies to specific bond lengths. Secondly, within the laser 
focal spot, there is a spatial variation of intensity which has to be taken into account above 
the saturation intensity. Thirdly, while present calculations only consider parallel electronic 
transitions, we must consider results averaged over molecular orientation. Previous work by 
Plummer and McCann 1121 found that DC ionization rates decrease sharply as the angle of 
orientation of the the molecular axis with the field increases. The orientation dependence of 
dynamic tunnelling ionization has yet to be established. These refinements are likely to be 
more important in the very high intensity regime / ^ 10^^ W cm^^ rather than the regime 
I ^ 10^"* W cm~^. We intend to undertake refinements of our model to simulate these effects 
and produce accurate estimates of ion yields and ion energy spectra. 
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